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ilar to previous estimates. We also explore the possibility of a large g'NL as the reason for 
the observed excess of the massive galaxy clusters. This scenario, (7nl > 2 x 10^, appears 
to be in more agreement with CMB and LSS limits for the non-Gaussianity parameters 
and could also provide an explanation for the overabundance of large voids in the early 
universe. 
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Introduction 

Many models of inflation predict an unobservable non-Gaussianity of the primordial cur- 
vature perturbation C^. Thus a detection of primordial non-Gaussianity at any level would 
rule out whole classes of models. Although observationally the primordial perturbations are 
Gaussian to a great accuracy, there is nevertheless still much room for non-Gaussianities, 
which are usually parameterized by the lowest-order non-linearity parameters /nl and ^nl • 
They can be defined by an expansion around a Gaussian perturbation (^g with 

C = Cg + ^/nl(C| - (C|)) + ^5nlC| + O(Cg') • (1) 

This expansion assumes a specific local form of non-Gaussianity. The non-linearity param- 
eters have previously been constrained by both CMB and LSS measurements. The best 
constraints for /nl come from WMAP 7-year results |] -10 < f^l"^ < 74. For compar- 
ision, many inflaton models predict |/nl| ^ 0(1)- For q'nl the limits are less strict, with 
WMAP 5-year results giving -5.6 x 10^ < ^nl < 6.4 x 10^, while halo bias and LSS 
yield Q —3,5 x 10^ < ^nl < 8.2 x 10^. For the latter limits one assumes that /nl ~ 0. 
It is noteworthy that the estimates from CMB and LSS for ^^nl are comparable. The 
expectation is that the Planck Surveyor Mission should be able to limit /nl ^ C'(5) Q. 

Recently, considerations of primordial non-Gaussianity have been extended to studies 
of high mass galaxy clusters. The interest has been triggered by the fact that at least 15 
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high-redshift (z > 1.0) galaxy clusters have been observed with masses measured to be 
about ~ IO^^Mq [||. These clusters have redshifts in the range 1.02 < z < 1.62 while 
the central values of their masses lie the range 0.57 x IO^^Mq < m < 1.0 x IO^^Mq 
@, |, H |, |10|, 0, 0, |ll. It has been argued that in the standard ACDM-cosmology 
with purely Gaussian initial perturbations, the probability for the existence of such very 
high-mass clusters is diminishingly small 1 14 1 ^ . 

The situation changes in the presence of primordial non-Gaussianity. Since /nl and g'NL 
are the coefficients of products of several Gaussian variables, they modify the behaviour of 
the PDF in the tails of the distribution, non-zero /nl increasing the probability of massive 
clusters and decreasing the probability for large voids, and (7nl increasing the probability 
of both. Thus if interpreted as a primordial feature, the anomalous abundance of large 
massive clusters might be evidence for departure from the Gaussianity of the primordial 
perturbation js], [l^ . 

The effect of non-Gaussianity on the abundance of massive clusters is encoded in the 
cluster mass function n(M, z, /nl, 9nl)- However, deriving n is a very involved calculation 
already in the Gaussian case. When non-Gaussianities are present, certain problems arise 
in the use of the cluster mass function, which, if not accounted for properly, result in large 
errors in the final estimates for the non-linearity parameters. Here we present a careful 
examination of the influence of non-Gaussianities on the cluster mass function, resolving 
the discrepancies that exist in the literature, and finding what we believe is the correct 
observational lower limit on /nl from the observed number of massive clusters. Moreover, 
we extend the previous analyses to higher order statistics and find the cluster limit also 
on gNL- We argue that a large q'nl is in more agreement with the CMB and LSS limits 
and point out that it could also provide an explanation for the observed overabundance of 



large voids in the early universe (see [17, 19| and references therein). 



The paper is organized as follows. In section |l| we introduce the formalism necessary 
to perform the quantitative analysis. In section || we then derive analytical estimates for 
the skewness and the kurtosis of the spectrum of density perturbations for given /nl and 
^NL- After that, in section ^ we estimate the value of /nl required to explain the number of 
heavy clusters using both analytical estimates and numerical calculations. We also discuss 
the subtleties leading to unreliable underestimates of /nl found in the literature. In section 
^ we then repeat this analysis assuming that the contribution of /nl is insignificant, and 
that the non-Gaussian primordial statistics is dominated by i^nl- In section ^ we summarize 
our results and point out the virtues of large (7nl in providing a possible explanation for 
the apparent overabundance of voids in the universe. 



1. Method 



The theoretical Gaussian mass function was first calculated by using a spherical collapse 
model and was later improved by generalizing this to ellipsoidal collapse. Even so, 

""^It is also argued in |^5| that a more conservative treatment of survey volume and measurement bias could 
reconcile these clusters with ACDM; however this conclusion is only arrived at for the clusters individually 
and not the ensemble as a whole. 
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Cluster Name Redshift M200 10^ Mq Mass Reference 



WARPS J1415. 1+3612 


1.02 


00+2 83 
"^•"-•"^-1.80 


floli 

i 


bF 1 -CL J 2341-51 19 


1.U3 


5.40+2.80 


m 


CLJ1415. 1+3612 


1.03 


3-40ll]i!] 


111 


XLSSJ022403. 9-041328 


1.05 




i 


SPT-CLJ0546-5345 


1.06 


10.01^:^1] 


i 


SPT-CLJ2342-5411 


1.08 


rv nrv+1 80 





RDCSJ0910+5422 


1.10 


r. r>o4-3 70 

6.281^:™ 


|1C 


RXJ1053.7+5735(West) 


1.14 


2-ootJ:?[] 


1 -1 n 

|12 


AljoaJUzzoUo.UU4obzz 


1 00 
i.zz 


1 in+0.60 


1^ 


RDCSJ1252. 92927 


1.23 


2-0010:1° 


|1C 


RXJ0849+4452 


1.26 


Q 7n+i-90 

"J- '"-1.90 


lie 


RXJ0848+4453 


1.27 


1.8011:1° 


|1C 


XMMUJ2235. 3+2557 


1.39 


'-'"-3.10 


|12 


XMMXCSJ2215.9-1738 


1.46 


4.iol?:|o 


12 


SXDF-XCLJ0218-0510 


1.62 


0.5710:1^ 


|1| 



Table 1: A list of the 15 observed high-mass and high-redshift galaxy clusters that we use in the 
present analysis. This list was first compiled in 



the theoretical predictions do not match the results of simulations, and the best estimates 
for nQ,{M^z) still come from semi-analytical fits to N-body simulations. 

The case for the non-Gaussian cluster mass function is even more complex. The N- 
body simulations have been mostly performed for Gaussian distributions, and thus no 
N-body formula for general non-Gaussian distribution exists. If deriving the accurate mass 
function was difficult for the Gaussian case, the non-Gaussianity of the distribution makes 
the task even more cumbersome. Thus it is no surprise that even the best estimates for the 
non-Gaussian mass function are only roughly valid, and do not always match the N-body 
simulations even in the Gaussian limit (/nl — ^ 0, g'NL — ^ 0). 

Since the analytical derivations have difficulty including non-ellipsoidal collapse (how- 
|2l|, pi |2|, pi), and thus the non- Gaussian mass functions cannot be trusted 



ever, see 



directly, the ratio of the non-Gaussian to Gaussian mass functions is often used to 
estimate the effect of non-Gaussianity, 

TZ (M, z, /nl, 9nl) = yrr 7 ^ ^ , (1-1) 

so that the non-Gaussian mass function is given by the Gaussian mass function (e.g. a 
fit to the Gaussian N-body simulation) multiplied by the ratio, nNcC-^, -z, /nl, Snl) = 
nQ{M, z) TZ. This expression has the correct behaviour in the non-Gaussian limit, and the 
erroneous behaviour of the non-Gaussian estimates is furthermore assumed to cancel in the 
ratio. There is evidence from N-body simulations that although apparently ad-hoc, these 
assumptions do have some validity |3| . 

To estimate the value of the non-Gaussianity parameters, two functions are needed 
for the use of formula ( |1 . 1| ) : the numerical Gaussian mass function and the analytical 
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non-Gaussian mass function. For the non-Gaussian function there are several different 
expressions in the hterature, with the most common being the Matarrese- Verde-Jimenez 
estimate (hereafter called the MVJ estimate) . Other expansions include the Edgeworth 
expansion ^ and its generalisations [p^, ^ 23 1. Analysis of the applicability of 



the various expansions can be found in ref. |30| who combine MVJ with the methods of 



|21|, H, |3|. For the ranges of redshifts and masses that we need to consider here, however, 
the Edgeworth expansion no longer converges. For that reason, we also use the MVJ 
expressions for both the dependence on /nl and ^nl H • 



We will follow the MVJ [g^ convention for defining the ratio TZ{M, z, /nLjS'nl) by 

3 'S'3(/nl) , .4 5'4(s'nl)\ (iSec dSs , . 1 / 1 (^ec ^S. 



K = exp I S-l^ + ei-FJ 1^-;^ + 1^-^ + } (1.2) 



where 6ec = \/0.75 x 1.686 is the critical density of ellipsoidal collapse, 5s = y^l — decSs/S 
and ^4 = y^l — (5^^5*4/12. 53(M, /nl), 'S'4(-^, 5nl) and aM all defined in section |2|, are the 
smoothed, skewness, kurtosis and variance of the nearly Gaussian density perturbations, 
respectively. 

The other popular definition for TZ uses the Edgeworth expansion, but is not suit- 
able for our purposes for reasons described in |^]. Specifically, the Edgeworth expansion 
involves an expansion over the terms appearing in the exponentials above. Despite the 
inherent smallness of 5*3 and ^4, the terms involving = Sec/crM can (and do) become 
large enough to overcome this. The Edgeworth expansion convention has been well tested 
in N-body simulations [p5| , |26| and performs well; however we will need to use our mass 
function beyond the ranges of these simulations. In [23| a full theoretical, non-Gaussian 



mass function is accurately determined that returns the Edgeworth expansion for 7^(/nl); 
however it is explicitly assumed that v^auS^ <C 1 to arrive at this result, making the result 



correct, but not applicable to the very high mass clusters. In |3C] this condition is relaxed 
and the full non-Gaussian mass function is calculated taking this quantity into account 
non-perturbatively. The TZ that results from this mass function is equivalent to eq. ( |1.2| ) 
up to the accuracy with which we know the masses of the clusters studied in this work. 
It is expected that this mass function will break down when vcfmSs — 1, a condition we 
need to be aware of when integrating to high masses, but that does not hold for any of 
the masses and redshifts of the clusters themselves for /nl ^ 1000. For vumSs > 1 (i.e. 
/nl > 1000 and > 5) there are currently no N-body simulations or theoretical results 
available. To derive upper bounds on the value of /nl allowed by clusters this problem will 
need to be addressed in the future. 

From ra(M, 2, /nl, 5nl) we then calculate the expected number of clusters of a given 
mass range, over a given redshift interval and with a given sky coverage using the following 
formula: 

rzf /.A/max 

Exp(Mi.ange,2;range,/NL,5'NL) = / dz dM fsky — n^Q, (1.3) 

m 

where /sky is the fractional sky coverage and (W{^z)ldz is the volume element at redshift z. 
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For both /nl and 5nl we quote results using the CMB convention. We also use 
WMAP5 g cosmological parameters^ : h = 0.705, = 0.96, = 0.28, Qbh^ = 0.0227 
and as = 0.812. 



2. Analytical estimate for Ss and ^4 

Due to a variety of differing claims in the literature, before performing the numerical 
analysis to extract information relating to /nl and ^nLi it will be useful to have some 
approximate, analytical, estimates for the effects both these parameters will have. We start 
by estimating the sizes of aS^ and These quantities are central to our calculation 

and different values have been claimed for them in the literature using identical definitions. 



2.1 Estimation of an and an 

Before we can estimate and 6*4 we will first estimate the variance of the density pertur- 
bations smoothed over a radius R. This will be useful for testing our approximations and 
will make the eventual estimations of 5*3 and 6*4 much more straightforward. The variance, 
(Tr, is defined by 

2 _ f°° 2 



aj,= I -ai{k,z)V{k), (2.1) 







where V is the almost scale invariant primordial power spectrum of and 

aR{k,z) = -^D{z)(^) T{k)WR{k). (2.2) 



Hq is the Hubble rate now, D{z) is the linear growth function, T{k) is the transfer function. 
For all results in this work we use the transfer function of [Bsl with the modified shape 



parameter of |34|. The window function is given by 



which is the Fourier transformation of a window function that is a top-hat in real space. The 
window function Wfj{k) is peaked at A; = 0; however the combination k'^WR{k) in eq. ( ^^ ) 
is closely peaked around kR ~ 1 with an amplitude tx 1/R^. T{k) is approximately one 
for the largest scales, but decreases for smaller scales that re-entered the horizon during 
the radiation dominated era. It does not decrease fast enough on these scales to stop the 
peak of the whole integral from still occurring at kR ~ 1. 

It is customary to define the normalisation of the power spectrum 'P{k) through the 
parameter a%, which is just evaluated at R = 8/i^^Mpc and 2; = 0. Given this and the 
above, it is possible to write, to a reasonable approximation, the following expression for 



^We do not use WMAP7 cosmology because the N-body simulations to date have only been performed 
using WMAP5 parameters. The deviations in our results arising from this will be much less than the errors 
implicit in the method itself. 
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an = aR,JV{l/R)D{z) (2.4) 

= J'-X^4^Diz), (2.5) 
^ \R) T(1/8) ^ ^ ' 

where is defined as aij(l/ii, 0). Ostensibly eq. ( |2.4D defines cr/j as a function of a/j and 
V] however, when estimating 53 and £'4, we will find it more useful to use this equation 
as a means to replace ur by the two quantities ctr and V that we know the approximate 
sizes of. 

In figure |^ we plot eq. ( p.5p as well as a numerically calculated curve for a range of 
R values at 2; = 0. Both curves will change identically with redshift through the linear 
growth function D{z). 

2.2 Estimation of crS^ and a^S^ 

Using the same arguments as above it is possible to estimate expressions for and 54. 
We start with the skewness which is given by the following expression, using /nl defined 
in eq. (|l|), 

a%S3{R) = 3hL ^aRiki)Vi -^aR{k2)V2 d^LaR{k^) , (2.6) 
Jo ki Jq k2 y_i 

where k"^ = k\ + k2+ 2jjLkik2 and Vi is short for V{ki). 

The three factors of aR{k) above will scale in the same way as they do for ctr. That 
is they will be highly peaked around the value k = 1 /R, with magnitudes proportional to 
The integrals over ki and k2 can be evaluated in exactly the same manner as the k 
integral for aR. The integral can be evaluated by appealing to the fact that ki and k2 will 
be ~ 1/R but not always exactly equal to 1/R. Therefore, over the full range n = —1 to 
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1, it will always be possible to find values of ki ~ 1/R and k2 — 1/R such that ^3 ~ 1/R. 
Therefore we can proceed by also substituting aR{k^) with an and multiplying the rest of 
the expression by 2 (the range of the 11 integral). This gives, 

4^3 (i?) = Qhi^alV^ = 6/nl4^'/'. 

Therefore, to the degree of our approximation, the quantity aiiS^{R) is scale independent 
and given by a very simple formula. For = 0.812, V ~ 2.4 x 10^^ which allows us to 
make the following estimate, 

anSsiR) ~ 3 x IO'^nl- (2.7) 

The growth factor D{z) is independent of scale, therefore the combination ajiS3{R) is 
exactly independent of redshift. 

A very similar argument applies for the kurtosis, Si{R). The kurtosis is given by,^ 



<7%S,{R) = -SNL n / T^««(^0^(^: 



\i=l ■ 

"1 rl f2n 

X I dfii dij.2 / dcpaRiki), 
-1 J-i Jo 



k\ = k\ + k\ + k1 + 2kik2fii + 2k2k^^2 



with, 



+2kiki ycos 0y 1 - /.if + /ii/i2j • 
If we follow the same process as we did for the skewness, this reduces to, 

a|S4(i?) = 245nl4^' = 245nl4^. (2.9) 
After we substitute the same value for V this gives, 

alSi{R) = 5.8 X lO^s^NL. (2.10) 

As with aiiS-i{R), this is also exactly independent of redshift. 

It was stated in an earlier version of this paper that our values of aRS^^R) and a'j^S4{R) 
were similar to those given in refs. but very different to the ones given in ref. [p9| . 



The authors of ref. |2£] have revised their paper and all sets of results are now in agreement. 
The value of 5*3 used in H] also matches our numerical result and analytic estimate.'* In 
our later numerical calculations we calculate Ss{R) fully numerically. We do not do the 



same for S4{R). The implications of this are discussed briefly in Section 4.2. Our numerical 



calculation of S3 matches our analytic estimate to within a factor ~ 1 — 1.5. 



^This particular expression for the kurtosis was first written down in an earlier version of ref. | 
*Ben Hoyle, private correspondence 
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3. Estimating /nl 



It was asked in ||5|, |T6| "what is the probabihty that a given cluster is the 'most massive' 
cluster in the survey window?". Here, the phrase 'most massive' strictly speaking means 
least probable because the redshift dependence of the mass function is always included. 
That is, a less massive cluster that collapsed much earlier could still be the 'most massive' 
cluster by this definition. To answer this question both references create three mass bins. 
One mass bin contains the central mass of the cluster as well as all the masses contained 
within the 1-a error range either side of it. The second mass bin contains all the masses 
above the upper bound of this previous range. The probabilities are obtained by first 
calculating the expected number of clusters observed in each bin, then Poisson sampling 
from each of these distributions a large number of times (10^) and asking which bin contains 
the largest cluster each time. It is assumed that the third bin, containing masses beneath 
the 1-a error range always contains at least one cluster in any given survey. 

Both refs. |5|, saw a sudden change in their calculated probabilities at a particular 
value of /nl- This point occurs at different values of /nl in each reference. The main 
symptom of this change is a sudden leap in the probability that the upper mass range 
contains at least one observed cluster. The references interpret this as being caused by the 
mass around the lower limit of this range becoming more probable at this value of /nl- 
It is striking that this transition occurs so rapidly in each reference and at considerably 
different values of /nl- We believe the true cause of this effect is the breakdown of the 
mass function being used in each reference at a much larger mass, corresponding to the 



arbitrary upper bound in their numerical calculation of the integral in eq. (1.3). 



We numerically examine this effect in section |3.2| ; however in section |3.1| below, fol- 
lowing the spirit of section ^, we attempt to give an approximate analytic estimate of how 
large /nl would need to be to observe the magnitude of effect seen in refs. [|, |16|. 



3.1 Estimate of /nl required to see the effect in ref. 

The cluster XMMUJ2235. 3+2557, with mass 7.7tti x IO^^Mq, observed at redshift 
z = 1.39 (see table ||) is found by |^ and us to be one of the two least probable clusters. 
This cluster was detected by an X-ray survey |35]. The full survey footprint of various 



X-ray surveys is estimated in M to be 283 sq. degrees. If we use the Jenkins et al mass 



function of ref. 1 36] in eq. (p..3|), and then integrate from redshift z = 1.39 — )• 2.2 and from 
mass M = 12 X 10^^ Mq — )■ oo we find the expected number of clusters in this upper bin 
should be 0.0019. This translates into a probability for this mass bin of ~ 0.002 which is 
consistent with figure 2 of ref. j^. 

/nl will change this probability through the ratio in eq. (|1.2D which to a good approx- 
imation is, 

1.3 



n = expi^-{aSs)j (3.1) 

where we have made the useful substitution, = bed^M to simplify the result. To the 
degree of approximation we are going to make, all the scale and redshift dependence is 



found in v and aS'^ is given by eq. ( |2.7| ) 
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We need to estimate the value of at M = 12 x IO^^Mq and redshift z = 1.39. To do 
this, we will use eq. (|2.5| ) for cjm (note that 6ec — 1.46). Eq. ( |2.5D gives aji, therefore we 
need to find what smoothing scale R{M) corresponds to each mass. We are interested in 
the comoving smoothing scale, therefore R{M) is given by 

f 3M , , 

R = , (3.2) 



where pm is the present density of matter in the Universe. For a critical density of pc = 
2.775//1 X 1O"M0 (/iMpc-^)3, h = 0.7 and Q.m = 0.28, 

R ~ 6M|f (/j-^Mpc), (3.3) 

where M14 is the mass of the cluster in fractions of IO^^Mq. Therefore, for M14 = 12, we 
arrive at i? ~ 13.7. Upon direct substitution into eq. (|2.5[) this gives for 



z.(Mi4 = 12,z) = -^. (3.4) 

Finally, one obtains for D(1.39) = 0.53 that u = 5.3. We find v = 5.0 using our numerical 
code. 

We are seeking to find the value of /nl necessary to increase the probability of this mass 
bin to a value indistinguishable from one. To be conservative we will stop at P = 0.95. For 
a Poissonian distribution to have a probability of 0.95 that there is at least one event, the 
expected number of events will need to exceed 3. Thus, we need an /nl that will increase 
the expected number of clusters belonging to the mass bin by a factor of 3/0.002=1500. 

If we put this into eq. (^), we find 

ln(1500) = ^3 X 10-Vnl, (3.5) 

D 

which gives /nl — 1000, a factor of two larger than that found in ref. Q. Although finding 



a value of /nl much closer to this, ref. [16| are not safe. This is due to them having a much 
smaller survey window of 11 sq. degrees. This will result in a much smaller probability 
for the cluster to exist in the Gaussian case, hence increasing the value of /nl required to 
make this mass bin probable. 

3.2 Discussion of Gaussian mass functions 

Following ^ 16], we have used above the Gaussian mass function from Jenkins et al. p6| . 
The Jenkins et al. result used in ref. ||5| is a fit to N-body simulations, which corresponds 
to ACDM and the spherical overdensity group finder, and is given by 

no{M,z) = lf(-''-^\ (3.6) 



Af V dlnM 

with 



/ = 0.301 exp [-1 \n(Tll{z) + QMf^"^] . (3.7) 
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In ref. |3£] this (and other) mass functions were tested and found to match well to 
simulations up to a value of Ino""^ < 1. At redshift z = this corresponds to a mass of 
5 X 10^^ Mq/H, which is well above the mass of any clusters we are considering. However, 
at redshift z = 1.4, a mass of IO^^Mq would correspond to In cr"^ = 1.2. As one probes 
deeper into redshift space, this mass function needs to be used at values of lno"~^ further 
and further beyond those for which it is tested. 

This does not necessarily mean the mass function will be unusable outside of this range. 
The clusters themselves are not far outside of the tested range of these mass functions so 



it could be assumed that any errors introduced will be small if one applies eq. (1.3) up 



to these masses and redshifts, but no further. However the redshift range probed by the 



cluster surveys is assumed in both refs. 16] to extend to z = 2.2. Moreover, when 



considering the upper mass bin described in section 3T, in principle one is required to 
integrate the mass to infinity. In practice this means integrating to some high mass where 
the full cluster mass function is assumed to be negligible. For this assumption to hold, it 
is imperative that at these redshifts and masses the Gaussian mass function decreases fast 
enough with increasing (7~^, since the non-Gaussian contribution is growing exponentially, 
with 7^ ~ exp(10~^/NLO-"^). 

Once either /nl or a^^ grows large enough eq. (3^) will not decrease fast enough. 
In figure ^ we have plotted the full mass function at redshift z = 2.2 for three values of 
/nl- The mass function is weighted by the volume element dV/dz to give the values on 
the y-axis more intuitive meaning. Each curve is effectively a plot of the integrand of 
eq. ( |1.3| ) against mass for various /nl values. It is clear that the formalism breaks down 
at a particular mass for a given redshift and /nl- The result of this breakdown is a very 
rapid turn-up in the full mass function at this mass. It is this effect and not a physical 
increase in the probability of clusters themselves existing that causes the sudden turn-up 
in the figures of refs. I, H]. 



In ref. Q, the upper limit of the integral in eq. (L3) was set to 1.5 x 10^^Mq//i.^ They 
also see the sudden turn-up in probability at /nl — 500. On figure |2| we have drawn a line 
upwards at M ~ 1.5 x 10^^ Mq/Ii. It is very clear from where it crosses the /nl = 500 curve 
that /nl = 500 is the first /nl value where this breakdown in the mass function will be 
seen at this upper limit. The different value of /nl in ref. [^6| for seeing this turn-up will be 
a result of a different upper limit used in eq. (|1.3|). For context, a mass of 1.5 x IO^^Mq/Zi 
at a redshift of 2.2 corresponds to a value of a^^ = 10.8, which is well outside the range 
quoted as valid for the Jenkins mass function, eq. (p.?]). 

To avoid this problem it is necessary to either use a Gaussian mass function that scales 
correctly at these high masses and redshifts, to use a full non-Gaussian mass function in- 
stead of the ratio method of eq. ( |l.lD or to cut off the integral before this breakdown occurs. 
Unfortunately no N-body fits have been tested this far into the tail of the distribution in 
both mass and redshift. In the absence of a tested function, we can ask what asymptotic 
behaviour we expect from theory. Thankfully, although differing in exact functional form, 
all theoretical models reproduce the same behaviour for the Gaussian mass function at large 



^Private correspondence. 
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Figure 2: Evidence of tlie brealcdown of tlie non-Gaussian mass function at large masses wlien 
using the ratio method and the Jenkins et al. Gaussian mass function. Plot is of dV/dz x n(M, z) 
versus mass, at redshift z — 2.2 and for some characteristic values of /nl- 



masses and redshifts, with / ~ exp(— c/cr^) in equation (3.7). The form of eq. ( |3.7D makes 
it very difficuh to judge how it will behave asymptotically as compared with exp(— c/cx^). 



Therefore, we will use the formula from Tinker et al. |37| , given by 



+ 1 



e""/"' , (3.8) 



where for an spherical overdensity A = 200 they find A = 0.186, a = 1.47, h = 2.57, and 
c ~ 1.19. This numerical fit also matches the most sophisticated theoretical models over a 
range of = 0.5 — t- 3, to within the errors introduced in the theoretical models [p2| . 

Figure ^ is a comparison of both mass functions. At smaller masses they match closely 
but diverge from each other at larger masses. This indicates that the Jenkins et al. mass 
function cannot be behaving asymptotically as / ~ exp(— c/o"'^). This is the origin of the 
problematic behaviour seen in figure |2[ In figure Q we plot the analogue of figure ^ using 
the Tinker et al. mass function. Two things are immediately apparent. Firstly, close to the 
relevant cluster masses it is much better behaved than the full mass function derived using 
the Jenkins et al. Gaussian — this is exactly what we expected from figure |3|. Secondly, it 
also breaks down if probed to large enough masses. This indicates that even a mass function 
that asymptotically matches the best theoretical models at both the Gaussian and non- 
Gaussian limits, breaks down at large enough masses. This too, however, is not unexpected 
and will occur even for the full theoretical non-Gaussian mass functions that do not use 
the ratio method of eq. (|1 . 1| ) . This is discussed briefly in [^]. This breakdown occurs at 
the point that 1^(0" 5*3) ~ 1 which is the point where the non-Gaussian contribution to the 
density perturbations becomes comparable in size to the Gaussian contribution. For the 



^This very simple formula matches the more complicated red-shift dependent formula presented in the 
update Ipsl] to within 5%. We use this formula due to its clear and desired asymptotic behaviour for small 
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Figure 3: Comparison of the Gaussian 
mass functions from Jenkins et al. [ p6[ and 
Tinker et al. [|3 plotted against mass, for 
z = 1. The lower plot is simply a zoomed in 
version of the upper plot. 
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Figure 4: Evidence of the breakdown of 
the non-Gaussian mass function at very large 
masses using Tinker. Plot is of dV/dz x 
n{M, z) versus mass, at redshift z = 2.2 and 
for the same /nl values as fig. ^. 



least probable clusters so far detected, /nl — 1000 is sufficiently large for the breakdown 
to occur. To set an upper bound on /nl in the future it will be necessary to extend the 
non-Gaussian mass function to v{aS-i) > 1. 

With a clear understanding of the limitations of the two mass functions, it is pertinent 
now to perform the analysis of section numerically. The results are presented in figures 
III, |6| and |^ for three of the clusters in table |l|. Figure |^ is the correct plot using the Tinker 
et al. mass function. By comparing figures I and I it is clear that, with the Jenkins et 
al. mass function, it is possible to cause the sudden, unphysical, jump in the probability of 
the upper mass bin to occur at any arbitrary value of /nl- This would be done by tuning 
the cutoff of the integral in eq. ( |1.3D to cause the unphysical turn-up of figure ^ to creep 
into the integrated range at precisely the chosen value of /nl- More positively however, it 
is also clear from examining all three figures, that until this unphysical jump does enter 
the integral, the results are consistent. This gives us confidence that the results presented 
in refs. ||, 16 1 will be correct whenever they correspond to /nl beneath this critical of /nl 
for their respective choice of Mmax- For rsf. [^J this will be /nl ^ 500, for ref. [16[ this will 
be /nl < 700. 

3.3 Robust numerical results for /nl 

We now seek to use the entire ensemble of 15 clusters to constrain /nl- For results presented 
in this section we always calculate the skewness, 5*3, fully numerically. We follow the 
method of ref. very closely including using their conservative estimates for /sky- For 
/nl < 500 we find similar results. We denote Mj as the random variable giving the mass of 
cluster i and we denote M,, af and as the mass, upper error and lower error respectively 
of cluster i quoted in table |l|. We take In Mj to be normally distributed with a mean of 

^ii = \nMi. (3.9) 

Also, we take the standard deviation, cjj, of InMj to be 



2 \Mi- aT 



(3.10) 
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Figure 5: The probability that each of the three labelled clusters could exist and be the "most 
massive" cluster in the survey using the Jenkins et al. mass function. Following the convention 
of the solid line depicts the probability that the "most massive" cluster in the survey is more 
massive than the labelled cluster, the dotted line depicts the probability that the "most massive" 
cluster's mass falls within the 1-a error range of the labelled cluster and the dashed line depicts the 
probability that the "most massive" cluster is less massive than the labelled cluster. For this curve 
we set Minax = 1.5 X 10^^ Mq/H in equation (O). 
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Figure 6: This is an identical figure to figure ||, except we set Af,„ax = 4.5 x lO^^MQ/h in eq. (1.3) 
This is clear evidence of the cutoff dependence of the results due to the effect seen in figure |[ 



For each of the 15 clusters, we then sample Mj 10^ times from this distribution. For each 
of the 10^ mass samples we calculate the total number of clusters expected in the survey 
window at equal or higher mass and redshift. We then Poisson sample over each one of these 
expectation values. The total number of these samples that return a value greater than zero 
allows us to form a probability, Pi, that the given cluster can exist, marginalised over its 
assumed log Gaussian distribution. We then multiply the individual cluster probabilities 
together to form the probability, P{f^L) = Y\Pij that the full ensemble could exist in 
the survey window. Finally, we repeat the above analysis for increasing values of /nl and 
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Figure 7: This figure is identical to figures || and g except it has been calculated using the Tinker 
et al. mass function. These results are insensitive to the choice of Mmax- 



1 




Figure 8: This figure is an exact analogue of figure ^ but for (/nl- The solid curves depict the 
probability that the 'most massive cluster' is greater than the labelled cluster, the dotted curves 
depict the probability that the 'most massive' cluster fits within the 1-a error range of the labelled 
cluster and the dashed line depicts the probability that the 'most massive' cluster is less massive 
than the labelled cluster. The appropriate labels are present on figure ||. 

record the value of /nl when -P(/nl) = 0.05. 

In figure ^ we have plotted ^(/nl)- If we compare this to figure 4 of ref. Q we see 
that the two results match well until /nl — 500. From the data used to produce this figure 
we find a lower bound on /nl|p=o.05 with WMAP 5 year parameters of /nl|p=o.05 > 412 
required to make the existence of this ensemble of clusters consistent with ACDM. ref. Q 
quote the corresponding value to be /nl|p=o.05 > 476. This is close, but perhaps more 
different than we would initially expect. We suggest the differences will come from different 
definitions of ^ and a; although it could also potentially come from our use of ref. [^]'s 
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Figure 9: The probability that the ensem- 
ble of clusters in table ^ could exist as a func- 
tion of /nl- 
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Figure 10: The probability that the en- 
semble of clusters in table ^ could exist as a 
function of ijnLj with /nl ^ 50. 



fitting formula for the transfer function, compared to ref. who numerically integrate 
theirs using the icosmo package [^]. Ordinarily one would not expect > 10% errors; 
however a systematic < 10% error in all 15 clusters could accumulate into a discrepancy 
of this size in -P(/nl)- 

Our main motivation is to compare this lower bound for /nl to the one we will derive 
in section ^ for (^nl- Therefore, for concision we will not marginalise this /nl constraint 
over the cosmological parameters. The results would be similar to that obtained by ref. ||5|, 
except for regions where /nl > 500, as discussed earlier. 

The analysis in ref. [16| uses a slightly different method to calculate their final con- 
straints on /nl, but gain a similar value, quoting /nl = 449 it 286. Their quoted numbers, 
however, are more dependent on the spurious behaviour caused by the breakdown of the 
mass functions, so their results change more significantly when this is removed. Never- 
theless, the effect of removing it would be to make larger /nl more probable. This would 
increase the central value of /nl that they quote, making their conclusions about non-zero 
and scale dependent /nl even stronger. Both a marginalisation over cosmology and a 
re-analysis similar to ref. |l^ would be interesting future calculations. 



4. Analytic estimate and numerical results for gfNL 

We now repeat the calculations of the previous section, but for (7nLi assuming /nl is small 
enough to be negligible in comparison (in practice this means /nl ^ 50). A positive /nl 
would mildly reduce our quoted lower bound on (^nl and a negative /nl would mildly 
increase it. 

4.1 Estimate of ^nl required to give similar effect as /nl 



To estimate the effect of (7nl we use the (7nl dependent part of eq. (1.2), which in terms of 
V, is to a good approximation 

7^:^expf^(a254)y (4.1) 
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This equation appears in the calculation for (/nl in exactly the same manner as eq. 
does for /nl- Therefore, in order for a given (/nl to cause a similar sized effect to a given 
/nl it is necessary that 

^(^^54) '^{aSs) (4.2) 

or, 

2 X 10^ ^ , ^ , 

ffNL = /NL- (4.3) 

For our least probable cluster, with u ~ 4.5, this gives (7nl — 4.4 x 10^ /nl- The least 
probable clusters dominate the departures from ACDM, therefore to a reasonable approx- 
imation we can use this result to relate our earlier constraints on /nl to (^nl- We found 
earlier that /nl — 410 was sufficient to give the ensemble of clusters a probability of 0.05 
of existing. We therefore expect for ^nl — 1-8 x 10^ the same to be true. As we will see, 
this is indeed the case. 

4.2 Numerical qtssl, 

Here we present the numerical results for (/nl using precisely the same methodology that 



was used for /nl in section 3^. In ||3| it was found that in the mass range 10^^ ^ M < 
5 X lO"*^^ Mq//i, varies in the narrow range ~ 4 — 6 x 10~^gNL (note that this particular 
combination, a'^S^, is independent of redshift). This matches closely to our analytical 
estimate, a'^S^ ~ 5.8 x IO^^^nl- Therefore, we choose to trust the results of ref. and 
take CT^54 = 5 X IO'^^nl to be constant over all redshifts and the mass range 10^^ ^ M < 
5 X 10^^ Mq/H. This range comfortably encompasses all of the clusters in table |l|. We 
expect this approximate method of defining (7nl should introduce errors in (7nl of < 20%. 

Due to the more rapid scaling of TZ with for ^4 as compared to ^3, we noticed that 
at z = 2.2 and with M^ax > W^^Mq, the full mass function sometimes broke down within 
the range of our integral. This occurs even with the Gaussian mass function from Tinker 
et al. [^]. This is due to ^'^{a'^S^) becoming of order one, rather than due to a breakdown 
in the ratio method. To avoid this affecting our results, we imposed a cutoff in our integral 
at v = 7.2. This is comfortably larger than values near any of the cluster masses and 
redshifts. It is also small enough that the mass function does not break down for any 
interesting values of q'nl- We tested the robustness of this approximation by varying the 
arbitrary cutoff. Within the ranges quoted and plotted, our results did not change. The 
changes that did occur were for very large (7nl precisely where we would expect the mass 
function to break down near u = 7.2. 

In figure ^ we show the plot of ^nl that is equivalent to figure 0. We see very similar 
behaviour to the plot of /nl- We would also see similar behaviour to figures |5| and ^ 
if we were to use the Gaussian mass function from Jenkins et al. In figure |lO| we show 
the plot of (7nl that is equivalent to figure |^. From this figure we see that our estimate of 
ffNL — 1.8 X 10^ is very accurate. When we extract the precise value at which -P(gNL) = 0.05, 
we find g'NL|p=o.05 > 2.0 x 10®. If we compare this lower constraint to the upper constraints 
listed in the introduction, we note a small degree of tension. However the amount of tension 
is much smaller than in the /nl case. 
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Figure 11: The effect of varying a% on our ^nl estimate (i.e. the value of ^nl for which the 
ensemble of clusters in table ^ exists with probability P — 0.05). 



Finally, in figure 11, we give a plot of gNL|p=o.05 against cxg. To calculate the effect 
of changes in as on ^4 we multiplied by the ratio (t|/0.81^ (note: erg = 0.81 was the 
value for which ref. calculated their value for cr^S'4). This takes into account the effect of 
(Tg on 5*4 exactly. This is because cjg only defines the normalisation of V and does not enter 
into any scale dependent quantities in the integral defining 54. The important result of 
this figure is that for large enough erg, but still within observational constraints, 5nl|p=o.05 
itself can be brought within current observational constraints. This would suggest that we 
can fully explain the existence of these massive, high-redshift clusters without appealing 
to scale dependent non-Gaussianity. This is not the case for /nl- 

This 'suggestion' might in fact be true, but we caution against taking it too seriously. 
The constraints on (^nl quoted in the introduction were obtained for a value of (Tg ~ 0.81. 
If we increase erg, then the degeneracy we see here that allows a smaller g-^i^ will force a 
smaller (7nl in those constraints as well. The open question is which degeneracy is more 
sensitive. This would be an interesting topic to pursue in the future. 

5. Conclusions 

We have found and remedied a common problem in the literature |l6| that occurs 
when estimating /nl from massive high redshift clusters. We found that the value of 
/nl necessary to make the existence of these massive clusters certain (i.e. to give 100% 
confidence that they would be detected in the given surveys) is /nl ^ 1000, which is 
significantly larger than what has been claimed (e.g. /nl > 550 in [^). We have also 
demonstrated that this discrepancy is explained by the fact that |^ and |jl6[ use a non- 
Gaussian mass function, derived from a Gaussian mass function that is not valid in the 
whole range of In cr~^ needed for the computation. Most importantly the Gaussian mass 
function does not scale correctly outside of this range. We have rectified this oversight 



by using a Gaussian mass function by Tinker et al. [37] that has the proper asymptotic 



behaviour. We also found that below a certain critical /nl value, the results in both refs. ||5|, 
can be trusted with confidence. However, this value differs for each references and 

depends on an arbitrary cutoff of the mass integral in equation ( |1.3D . Our 95% confidence 

lower bound of /nl ^ 410 agrees reasonably well with reference Q because it occurs below 

this critical /nl value for their cutoff. 

We then considered the case where the non-Gaussian part of the primordial spectrum 

is dominated by i^nl- Such a situation can easily arise e.g. in curvaton models p2| . 

We estimated (7nl > 2.0 x 10^. We thus demonstrated, that within current observational 

limits, fl'NL appears to have more potential to explain the observed excess of high-redshift 

massive clusters. 

Non-Gaussianity is not the only potential explanation for the tension caused by the 
existence of these high redshift clusters. A systematic error that consistently over-estimated 
the masses would have a similar effect. In |5| this possibility was considered. To minimise 
the effects of systematic errors, the mass measurements used were always taken to be 
those whose quoted errors were consistent with the smallest mass value. We use the 
same mass estimates. For many clusters, this involved comparing mass estimates from SZ 
effect measurements. X-ray measurements and weak lensing measurements. Any systematic 
errors would need to be present in all three measurement methods. It was also found in 
fsl that all the masses would need to be systematically over-estimated by 1.5 a in each 
measurement technique to make this ensemble of clusters fully consistent with /nl = 0. 

A different expansion history is another potential explanation for the existence of 
these clusters p3| , 44, |4^. A modified equation of state for dark energy that resulted in 
an earlier onset of the accelerated expansion would suppress structure growth at smaller 
redshifts/later times. This would have the effect of causing the linear growth function -D(z) 
to drop more slowly from low redshifts to high redshifts. The net result is more structure 
and thus more, large mass, clusters at high redshifts than what would be expected in 
ACDM. 

The estimates of /nl quoted here and in and ||l^ are far outside the observational 
limits set by WMAP. In Q and they suggested remedying this problem by introducing 
running of /nl- While this would explain the apparent discrepancy of the magnitude of 
/nl over different scales, it is also likely to introduce problems with the actual spectral 
index of the perturbations, since usually if /nl acquires running, so does the magnitude 



of the perturbations. (For discussions on running non-linearity parameters, see |4f:, 47, 48 



4|, ||.) 

We have demonstrated that instead of introducing non-zero n^^^^ the abundance of 
the massive clusters can be explained by introducing ^nl almost within the current obser- 
vational bounds. We also like to draw attention to a potentially important observational 
result: there appears to be an overabundance of large voids (see |17, 18| and references 
therein), and while large (7nl makes both heavy clusters and large voids more probable, a 
positive /nl actually makes the voids less likely, increasing the tension with observations 
even more [19|. In general it seems that a large value of (7nl is in much better agreement 
with all observations (CMB, clusters, halo bias, voids), than large values of /nl- With 
better measurements and N-body simulations a slightly smaller value of ^nl could perhaps 
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accommodate the abundance of the heavy clusters. This stresses the importance of ex- 
tracting hmits from the Planck microwave temperature map not only on /nl but also on 
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